#include "numb.h"
using namespace std;

/*
 * @brief Barrett模约减流程，根据相关预计算量计算  a mod q
 *
 * @param a 要约减的数a，长度为128bit
 * @param q 模数q,长度64bit
 * @param muHigh Barrett约减中的预计算量mu的高64位，mu为2^128−1除以n的商
 * @param muLow Barrett约减中的预计算量mu的低64位
 * @return uint64_t a mod q
 */
uint64_t SingleBarrettNew(__uint128_t &a, uint64_t q, uint64_t muHigh, uint64_t muLow)
{
    return Barrett128(a, q, muHigh, muLow);
}

/* *
 * 功能描述：Barrett计算流程的汇编优化，供上层Barrett函数调用
 * @__uint128_t& a:要约减的数uint128bit的数a
 * @uint64_t q：模数q
 * @uint64_t mu_high：Barrett约减中的预计算量mu的高64位，mu为2^128−1除以n的商
 * @uint64_t mu_low：Barrett约减中的预计算量mu的低64位
 *
 */
uint64_t Barrett128(__uint128_t &a, uint64_t q, uint64_t mu_high, uint64_t mu_low)
{
    uint64_t a_high = a >> 64;
    uint64_t a_low = static_cast<uint64_t>(a);
    uint64_t qq = 2 * q;
    asm volatile(
        // Step 1: Calculate a_low * mu_high, store high 64 bits in x1
        "umulh x1, %[a_low], %[mu_high]\n\t"

        // Step 2: Calculate a_high * mu_low and add to x1
        "umulh x0, %[a_high], %[mu_low]\n\t"
        "add x1, x1, x0\n\t"

        // Step 3: Calculate a_high * mu_high and add to x1
        "mul x0, %[a_high], %[mu_high]\n\t"
        "add x1, x1, x0\n\t"

        // Step 4: Calculate x1 * q
        "mul x1, x1, %[q]\n\t"

        // Step 5: Calculate a_low - x1
        "sub x1, %[a_low], x1\n\t"

        // Step 6: if (x1 >= qq) x1 -= qq
        "subs x0, x1, %[qq]\n\t"  // Calculate x1 - qq and set flags
        "csel x1, x0, x1, hs\n\t" // If x1 >= qq (hs), select x1 - qq, else
                                  // select x1

        // Step 7: if (x1 >= q) x1 -= q
        "subs x0, x1, %[q]\n\t"         // Calculate x1 - q and set flags
        "csel %[a_low], x0, x1, hs\n\t" // If x1 >= q (hs), select x1 - q, else
                                        // select x1

        // Store the result back to a_low
        // "mov %[a_low], x1\n\t"

        : [ a_low ] "+r"(a_low) // Output: a_low (modified)
        : [ a_high ] "r"(a_high), [ mu_low ] "r"(mu_low), [ mu_high ] "r"(mu_high), [ q ] "r"(q), [ qq ] "r"(qq)
        : "x0", "x1", "memory", "cc");

    return a_low;
}

/* *
 * 功能描述：对 64-bit 的 a 进行 Barrett 约减
 * @uint64_t a: 待约减的数
 * @uint64_t q：模数 q
 * @uint64_t mu：Barrett 约减的预计算量，mu = floor(2^64 / q)
 */
uint64_t Barrett64(uint64_t a, uint64_t q, uint64_t mu)
{
    asm volatile(
        // Step 1: Calculate a * mu, store high 64 bits in x1
        "umulh x1, %[a], %[mu]\n\t"

        // Step 2: Calculate x1 * q
        "mul x1, x1, %[q]\n\t"

        // Step 3: Calculate a - x1
        "sub x1, %[a], x1\n\t"

        // Step 4: if (x1 >= q) x1 -= q
        "subs x0, x1, %[q]\n\t"     // Calculate x1 - q and set flags
        "csel %[a], x0, x1, hs\n\t" // If x1 >= q (hs), select x1 - q, else
                                    // select x1

        : [ a ] "+r"(a)                // Output: a (modified)
        : [ q ] "r"(q), [ mu ] "r"(mu) // Inputs: q, mu
        : "x0", "x1", "cc"             // Clobbers: x0, x1, condition flags
    );
    return a;
}

/* *
 * 功能描述：对 64-bit 的 a 进行 Barrett 约减
 * @uint64x2_t a：存在NEON寄存器中的待约简的两个64bit数a
 * @uint64x2_t q：存在NEON寄存器中的待约简的两个模数q
 * @uint64x2_t mu：Barrett 约减的预计算量，mu
 * @uint64x2_t *res: 两个数的约简结果
 */
void Barrett64NEON(uint64x2_t *res, uint64x2_t a, uint64x2_t q, uint64x2_t mu)
{
    uint64_t res0 = 0, res1 = 0;

    asm volatile("umulh x1, %[a0], %[mu0]\n\t"
        "mul x1, x1, %[q0]\n\t"
        "sub x1, %[a0], x1\n\t"
        "subs x0, x1, %[q0]\n\t"
        "csel %[res0], x0, x1, hs\n\t"

        : [ res0 ] "=r"(res0)
        : [ a0 ] "r"(vgetq_lane_u64(a, 0)), [ q0 ] "r"(vgetq_lane_u64(q, 0)), [ mu0 ] "r"(vgetq_lane_u64(mu, 0))
        : "x0", "x1", "cc");

    asm volatile("umulh x1, %[a1], %[mu1]\n\t"
        "mul x1, x1, %[q1]\n\t"
        "sub x1, %[a1], x1\n\t"
        "subs x0, x1, %[q1]\n\t"
        "csel %[res1], x0, x1, hs\n\t"

        : [ res1 ] "=r"(res1)
        : [ a1 ] "r"(vgetq_lane_u64(a, 1)), [ q1 ] "r"(vgetq_lane_u64(q, 1)), [ mu1 ] "r"(vgetq_lane_u64(mu, 1))
        : "x0", "x1", "cc");

    *res = vsetq_lane_u64(res0, *res, 0);
    *res = vsetq_lane_u64(res1, *res, 1);
}

/* *
 * @uint64_t a：需要求相反数的值
 * @uint64_t& r：求相反数的结果
 */
void Negate(uint64_t &r, uint64_t a)
{
    r = -a;
}

/* *
 * @uint64_t& r:返回结果(a+b)%p
 * @uint64_t a:模加的操作数a
 * @uint64_t b：模加的操作数b
 * @uint64_t m：模数p
 */
void AddMod(uint64_t &r, uint64_t a, uint64_t b, uint64_t p)
{
    ComSubQ(a, p);
    ComSubQ(b, p);
    r = a + b;
    if (r > p) {
        r -= p;
    }
}

/*
 * @brief ：对两个64位整数进行加法并模p运算
 *
 * @param a：待加法操作的整数
 * @param b：加数
 * @param p：模数
 */
void AddModAndEqual(uint64_t &a, uint64_t b, uint64_t p)
{
    a = a + b;
    if (a > p) {
        a -= p;
    }
}

/* *
 * @uint64_t& r:返回结果(a - b)%p
 * @uint64_t a:模减的操作数a
 * @uint64_t b：模减的操作数b
 * @uint64_t m：模数p
 */
void SubMod(uint64_t &r, uint64_t a, uint64_t b, uint64_t m)
{
    if (m == 0) {
        return;
    }
    ComSubQ(a, m);
    ComSubQ(b, m);
    r = b % m;
    r = (a + m - r) % m;
}

void SubModAndEqual(uint64_t &a, uint64_t b, uint64_t q)
{
    uint64_t r = b % q;
    a = (a + q - r) % q;
}


/*
 * @uint64_t& r:返回结果(a * b)%p
 * @uint64_t a:模乘的操作数a
 * @uint64_t b：模乘的操作数b
 * @uint64_t m：模数p
 */
void MulMod(uint64_t &r, uint64_t a, uint64_t b, uint64_t m)
{
    if (m == 0) {
        return;
    }
    if (a > m)
        a = a % m;
    if (b > m)
        b = b % m;
    unsigned __int128 mul = static_cast<unsigned __int128>(a) * (b);
    mul %= static_cast<unsigned __int128>(m);
    r = static_cast<uint64_t>(mul);
}

/*
 * 功能描述：计算 a*b mod p，并且通过barrett约减来计算最终结果
 * @uint64_t& r：用来传递计算结果
 * @uint64_t a: 源操作数a
 * @uint64_t b: 源操作数b
 * @uint64_t p: 模数p
 * @uint64_t twop：模数p的两倍
 * @uint64_t prHigh: Barrett预计算量Pr的高64bit #在Context中定义对应的Pr矩阵
 * @uint64_t prLow:  Barrett预计算量Pr的低64bit
 */
void MulModBarrett(uint64_t &r, uint64_t a, uint64_t b, uint64_t p, uint64_t twop, uint64_t prHigh, uint64_t prLow)
{
    if (p == 0) {
        return;
    }
    if (a > p)
        a = a % p;
    if (b > p)
        b = b % p;
    asm volatile(
        // 计算 a*b  低位存x3 高位存x3寄存器
        "mul   x3 ,%[a] , %[b] \n\t "
        "umulh x4 ,%[a] , %[b] \n\t"
        // 计算Mul_l 与预计算量pr的乘积，高64位存在res中
        "umulh x1, x3, %[mu_high]\n\t"

        // 计算 a.high * mu.low + res :
        // 先计算高64位乘法，再累加res,要存放到res中
        "umulh x0, x4, %[mu_low]\n\t"
        "add x1, x1, x0\n\t"

        // 计算 a.high * mu.high + res,先计算低64位乘法，
        // 再累加res,要存放到res中
        "mul x0, x4, %[mu_high]\n\t"
        "add x1, x1, x0\n\t"

        // 计算 res * q 存回res
        "mul x1, x1, %[q]\n\t"

        // 计算 a.low - res,存回res
        "sub x1, x3, x1\n\t"

        // 计算res 和2*q的比较，如果直接减发生溢出就set flag
        "subs x0, x1, %[qq]\n\t"
        "csel x1, x0, x1, hs\n\t"

        // 计算res 和q的比较，如果直接减发生溢出就set flag
        "subs x0, x1, %[q]\n\t"
        "csel %[r], x0, x1, hs\n\t"
        // Store the result back to a_low
        //   "mov %[r], x1\n\t"
        : [ a ] "+r"(a), [ b ] "+r"(b), [ r ] "+r"(r)
        : [ mu_low ] "r"(prLow), [ mu_high ] "r"(prHigh), [ q ] "r"(p), [ qq ] "r"(twop)
        : "x0", "x1", "x3", "x4", "memory", "cc");
}
/*
 * 功能描述：计算 a 的barrett约减来计算最终结果
 * @uint64_t& r：用来传递计算结果
 * @uint64_t a: 源约简操作数a
 * @uint64_t p: 模数m
 * @uint64_t twop：模数p的两倍
 * @uint64_t prHigh: Barrett预计算量Pr的高64bit #在Context中定义对应的Pr矩阵
 * @uint64_t prLow:  Barrett预计算量Pr的低64bit
 */
void ModBarrett(uint64_t &r, uint64_t a, uint64_t p, uint64_t twop, uint64_t prHigh, uint64_t prLow)
{
    twop = 0;
    __uint128_t a_ = static_cast<__uint128_t>(a);
    uint64_t a_high = a_ >> 64;
    uint64_t a_low = static_cast<uint64_t>(a_);
    r = Barrett128(a_, p, prHigh, prLow);
}
void ModBarrett(uint64_t &r, unsigned __int128 a, uint64_t p, uint64_t twop, uint64_t prHigh, uint64_t prLow)
{
    twop = 0;
    uint64_t a_high = a >> 64;
    uint64_t a_low = static_cast<uint64_t>(a);
    r = Barrett128(a, p, prHigh, prLow);
}

/*
 * @uint64_t x：模逆的操作数x
 * @uint64_t m：模数p
 * @return  :返回结果x^-1%p
 */
uint64_t InvMod(uint64_t x, uint64_t p)
{
    int64_t m0 = p;
    int64_t x0 = 0, x1 = 1;
    int64_t q, t;

    if (p == 1)
        return 0; // 任何数 mod 1 都没有逆元

    while (x > 1) {
        // q 是商
        q = x / p;
        t = p;
        // m 是余数，a 是被除数
        p = x % p;
        x = t;
        // 更新 x0 和 x1
        t = x0;
        x0 = x1 - q * x0;
        x1 = t;
    }

    // 确保 x1 是正数
    if (x1 < 0)
        x1 += m0;

    if (x1 == -1) {
        //        throw invalid_argument("Inverse doesn't exist!!!");
    }

    return x1;
}
/*
 * @brief 计算x^y mod modulus
 * @param x
 * @param y
 * @param modulus
 * @return uint64_t
 */
uint64_t PowMod(uint64_t x, uint64_t y, uint64_t modulus)
{
    uint64_t res = 1;
    while (y > 0) {
        if (y & 1) {
            MulMod(res, res, x, modulus);
        }
        y = y >> 1;
        MulMod(x, x, x, modulus);
    }
    return res;
}

/* *
 * @brief ：计算 x 的逆元素（通过快速幂算法）。
 *
 * @param x：输入值 x。
 * @return uint64_t：返回 x 的逆元素。
 */
uint64_t Inv(uint64_t x)
{
    return Pow(x, static_cast<uint64_t>(-1));
}

/* *
 * @brief ：计算 x^y。
 *
 * @param x：底数 x。
 * @param y：指数 y。
 * @return uint64_t：返回 x^y 的计算结果。
 */
uint64_t Pow(uint64_t x, uint64_t y)
{
    uint64_t res = 1;
    while (y > 0) {
        if (y & 1) {
            res *= x;
        }
        y = y >> 1;
        x *= x;
    }
    return res;
}

/* *
 * @brief ：对输入值 x 进行位反转。
 *
 * @param x：输入值 x。
 * @return uint32_t：返回反转后的结果。
 */
uint32_t BitReverse(uint32_t x)
{
    x = (((x & 0xaaaaaaaa) >> 1) | ((x & 0x55555555) << 1));
    x = (((x & 0xcccccccc) >> 2) | ((x & 0x33333333) << 2));
    x = (((x & 0xf0f0f0f0) >> 4) | ((x & 0x0f0f0f0f) << 4));
    x = (((x & 0xff00ff00) >> 8) | ((x & 0x00ff00ff) << 8));
    return ((x >> 16) | (x << 16));
}

/* *
 * @brief ：计算 a 和 b 的最大公约数。
 *
 * @param a：输入值 a。
 * @param b：输入值 b。
 * @return uint64_t：返回 a 和 b 的最大公约数。
 */
uint64_t Gcd(uint64_t a, uint64_t b)
{
    if (a == 0) {
        return b;
    }
    return Gcd(b % a, a);
}
long Gcd(long a, long b)
{
    if (a == 0) {
        return b;
    }
    return Gcd(b % a, a);
}

/* *
 * @brief ：将给定的整数分解为其所有的素因数
 *
 * @param s：存储素因数的集合
 * @param number：待分解的整数
 */
void FindPrimeFactors(set<uint64_t> &s, uint64_t number)
{
    while (number % 2 == 0) {
        s.insert(2);
        number /= 2;
    }
    for (uint64_t i = 3; i < sqrt(number); i++) {
        while (number % i == 0) {
            s.insert(i);
            number /= i;
        }
    }
    if (number > 2) {
        s.insert(number);
    }
}

/* *
 * @brief ：查找给定模数的原根
 *
 * @param modulus：模数
 * @return uint64_t：模数的原根，如果不存在则返回 -1
 */
uint64_t FindPrimitiveRoot(uint64_t modulus)
{
    set<uint64_t> s;
    uint64_t phi = modulus - 1;
    FindPrimeFactors(s, phi);
    for (uint64_t r = 2; r <= phi; r++) {
        bool flag = false;
        for (auto it = s.begin(); it != s.end(); it++) {
            if (PowMod(r, phi / (*it), modulus) == 1) {
                flag = true;
                break;
            }
        }
        if (flag == false) {
            return r;
        }
    }
    return -1;
}

/* *
 * @brief ：查找模数下的 M 次单位根
 *
 * @param M：单位根的阶
 * @param mod：模数
 * @return uint64_t： M 次单位根的值，如果无法找到则返回 -1
 */
uint64_t FindMthRootOfUnity(uint64_t M, uint64_t mod)
{
    uint64_t res;
    res = FindPrimitiveRoot(mod);
    if (M == 0) {
        return 0;
    }
    if ((mod - 1) % M == 0) {
        uint64_t factor = (mod - 1) / M;
        res = PowMod(res, factor, mod);
        return res;
    } else {
        return -1;
    }
}

/* *
 * @brief ：测试一个数是否为素数
 *
 * @param p：待测试的数
 * @return bool：如果是素数返回 true，否则返回 false
 */
bool PrimeTest(uint64_t p)
{
    if (p < 2) {
        return false;
    }
    if (p != 2 && p % 2 == 0) {
        return false;
    }
    uint64_t s = p - 1;
    while (s % 2 == 0) {
        s /= 2;
    }
    for (long i = 0; i < 200; i++) {
        uint64_t temp1 = rand();
        temp1 = (temp1 << 32) | rand();
        temp1 = temp1 % (p - 1) + 1;
        uint64_t temp2 = s;
        uint64_t mod = PowMod(temp1, temp2, p);
        while (temp2 != p - 1 && mod != 1 && mod != p - 1) {
            MulMod(mod, mod, mod, p);
            temp2 *= 2;
        }
        if (mod != p - 1 && temp2 % 2 == 0) {
            return false;
        }
    }
    return true;
}

/* *
 * @brief ：计算给定整数 x 对 p 的 Shoup 表示
 *
 * @param x：整数
 * @param p：模数
 * @return uint64_t： x 对 p 的 Shoup 表示
 */
uint64_t x_Shoup(uint64_t x, uint64_t p)
{
    if (p == 0) {
        return 0;
    }
    unsigned __int128 temp = x;
    temp <<= 64;
    return uint64_t(temp / p);
}

/* *
 * @brief ：使用 Shoup 优化的乘法计算 a * b，结果模 q，并将结果存储在 res 中
 *
 * @param res：存储结果的变量
 * @param a：第一个数
 * @param b：第二个数
 * @param b_shoup：b 的 Shoup 表示
 * @param q：模数
 */
void MulModShoup(uint64_t &res, uint64_t a, uint64_t b, uint64_t bShoup, uint64_t q)
{
    unsigned __int128 v_w_mul = static_cast<unsigned __int128>(a) * bShoup;
    uint64_t hi = v_w_mul >> 64;
    res = a * b - hi * q;
    if (res >= q)
        res -= q;
}

/* *
 * @brief ：使用Barrett乘法算法对输入的多项式系数进行乘法并取模。
 *
 * 该函数通过NEON指令集实现Barrett乘法优化，执行大整数乘法并模一个指定的数
 * p。计算过程中使用了高效的硬件指令集，以加速模乘操作。
 *
 * @param r： 输出结果向量，存储计算后的模乘结果。
 * @param a： 输入的第一个多项式系数，表示待乘的第一个数。
 * @param b： 输入的第二个多项式系数，表示待乘的第二个数。
 * @param p： 用于取模的模数 p。
 * @param twop： 模数 p 的两倍值，用于进一步优化模操作。
 * @param prHigh： 预计算的模数 p 的高位部分，用于优化乘法运算。
 * @param prLow： 预计算的模数 p 的低位部分，用于优化乘法运算。
 */
void MulModBarrett_NEON(uint64x2_t *r, uint64x2_t a, uint64x2_t b, uint64x2_t p, uint64x2_t twop, uint64x2_t prHigh,
    uint64x2_t prLow)
{
    uint64_t res0, res1;
    asm volatile("mul x3, %[a0], %[b0]\n\t"
        "umulh x4, %[a0], %[b0]\n\t"
        "umulh x1, x3, %[prHigh0]\n\t"
        "umulh x0, x4, %[prLow0]\n\t"
        "add x1, x1, x0\n\t"
        "mul x0, x4, %[prHigh0]\n\t"
        "add x1, x1, x0\n\t"
        "mul x1, x1, %[p0]\n\t"
        "sub x1, x3, x1\n\t"
        "subs x0, x1, %[twop0]\n\t"
        "csel x1, x0, x1, hs\n\t"
        "subs x0, x1, %[p0]\n\t"
        "csel %[res0], x0, x1, hs\n\t"
        : [ res0 ] "=r"(res0)
        : [ a0 ] "r"(vgetq_lane_u64(a, 0)), [ b0 ] "r"(vgetq_lane_u64(b, 0)), [ p0 ] "r"(vgetq_lane_u64(p, 0)),
        [ twop0 ] "r"(vgetq_lane_u64(twop, 0)), [ prHigh0 ] "r"(vgetq_lane_u64(prHigh, 0)),
        [ prLow0 ] "r"(vgetq_lane_u64(prLow, 0))
        : "x0", "x1", "x3", "x4", "cc");
    asm volatile("mul x3, %[a1], %[b1]\n\t"
        "umulh x4, %[a1], %[b1]\n\t"
        "umulh x1, x3, %[prHigh1]\n\t"
        "umulh x0, x4, %[prLow1]\n\t"
        "add x1, x1, x0\n\t"
        "mul x0, x4, %[prHigh1]\n\t"
        "add x1, x1, x0\n\t"
        "mul x1, x1, %[p1]\n\t"
        "sub x1, x3, x1\n\t"
        "subs x0, x1, %[twop1]\n\t"
        "csel x1, x0, x1, hs\n\t"
        "subs x0, x1, %[p1]\n\t"
        "csel %[res1], x0, x1, hs\n\t"
        : [ res1 ] "=r"(res1)
        : [ a1 ] "r"(vgetq_lane_u64(a, 1)), [ b1 ] "r"(vgetq_lane_u64(b, 1)), [ p1 ] "r"(vgetq_lane_u64(p, 1)),
        [ twop1 ] "r"(vgetq_lane_u64(twop, 1)), [ prHigh1 ] "r"(vgetq_lane_u64(prHigh, 1)),
        [ prLow1 ] "r"(vgetq_lane_u64(prLow, 1))
        : "x0", "x1", "x3", "x4", "cc");
    *r = vsetq_lane_u64(res0, *r, 0);
    *r = vsetq_lane_u64(res1, *r, 1);
}

/*
 * @brief ：对两个64位整数进行加法运算，并对结果进行模运算
 *
 * @param r：存储结果的指针
 * @param a：第一个操作数
 * @param b：第二个操作数
 * @param p：模数
 */
void AddMod_NEON(uint64x2_t *r, uint64x2_t a, uint64x2_t b, uint64x2_t p)
{
    uint64x2_t sum = vaddq_u64(a, b);
    uint64x2_t cmp = vcgtq_u64(sum, p);
    *r = vsubq_u64(sum, vandq_u64(cmp, p));
}

/*
 * @brief ：对两个64位整数进行乘法运算，应用Shoup优化技术并对结果进行模运算
 *
 * @param r：存储结果的指针
 * @param a：第一个操作数
 * @param b：第二个操作数
 * @param bShoup：Shoup优化参数
 * @param q：模数
 */
void MulModShoup_NEON(uint64x2_t *r, uint64x2_t a, uint64x2_t b, uint64x2_t bShoup, uint64x2_t q)
{
    uint64_t res0, res1;
    asm volatile("mul    x3, %[a0], %[bShoup0] \n\t"
        "umulh  x4, %[a0], %[bShoup0] \n\t"
        "mul    x1, x4, %[q0]        \n\t"
        "mul    x0, %[a0], %[b0]     \n\t"
        "sub    x1, x0, x1           \n\t"
        "subs   x0, x1, %[q0]        \n\t"
        "csel   %[res0], x1, x0, lo  \n\t"
        : [ res0 ] "=r"(res0)
        : [ a0 ] "r"(vgetq_lane_u64(a, 0)), [ b0 ] "r"(vgetq_lane_u64(b, 0)),
        [ bShoup0 ] "r"(vgetq_lane_u64(bShoup, 0)), [ q0 ] "r"(vgetq_lane_u64(q, 0))
        : "x0", "x1", "x3", "x4", "cc");
    asm volatile("mul    x3, %[a1], %[bShoup1] \n\t"
        "umulh  x4, %[a1], %[bShoup1] \n\t"
        "mul    x1, x4, %[q1]        \n\t"
        "mul    x0, %[a1], %[b1]     \n\t"
        "sub    x1, x0, x1           \n\t"
        "subs   x0, x1, %[q1]        \n\t"
        "csel   %[res1], x1, x0, lo  \n\t"
        : [ res1 ] "=r"(res1)
        : [ a1 ] "r"(vgetq_lane_u64(a, 1)), [ b1 ] "r"(vgetq_lane_u64(b, 1)),
        [ bShoup1 ] "r"(vgetq_lane_u64(bShoup, 1)), [ q1 ] "r"(vgetq_lane_u64(q, 1))
        : "x0", "x1", "x3", "x4", "cc");
    *r = vsetq_lane_u64(res0, *r, 0);
    *r = vsetq_lane_u64(res1, *r, 1);
}

/*
 * @brief ：计算64位整数x的位数
 *
 * @param x：输入整数
 * @return 该整数的位数
 */
int BitsLen64(uint64_t x)
{
    if (x == 0)
        return 0;
    int len = 0;
    while (x != 0) {
        len++;
        x >>= 1;
    }
    return len;
}

/*
 * @brief ：根据Q和P的级别，计算RNS（残差数制）分解
 *
 * @param levelQ：Q的级别
 * @param levelP：P的级别
 * @return RNS分解后的级别
 */
int DecompRNS(int levelQ, int levelP)
{
    if (levelP == 0) {
        return levelQ;
    }
    return int(ceil(levelQ * 1.0 / levelP));
}

/*
 * @brief ：计算溢出余量，确保结果不会溢出
 *
 * @param level：级别
 * @param vec：待检查的数值数组
 * @return 溢出余量
 */
int OverflowMargin(int level, uint64_t *vec)
{
    __uint128_t pow64 = static_cast<__uint128_t>(1) << 64;
    uint64_t maxNum = vec[0];
    for (int i = 1; i < level; i++) {
        maxNum = max(maxNum, vec[i]);
    }
    if (maxNum == 0) {
        return 0;
    }
    return pow64 / maxNum;
}

/*
 * @brief ：基于2的幂分解，计算P的分解
 *
 * @param levelQ：Q的级别
 * @param levelP：P的级别
 * @param pow2Base：2的幂基数
 * @param qVec：Q的数值数组
 * @param pVec：P的数值数组
 * @return 分解后的P的级别
 */
int DecompPw2(int levelQ, int levelP, int pow2Base, vector<uint64_t>& qVec, vector<uint64_t>& pVec)
{
    if (pow2Base == 0) {
        return 1;
    }
    return (MaxBit(levelQ, levelP, qVec, pVec) + pow2Base - 1) / pow2Base;
}

/*
 * @brief ：计算Q和P的最大位数
 *
 * @param levelQ：Q的级别
 * @param levelP：P的级别
 * @param qVec：Q的数值数组
 * @param pVec：P的数值数组
 * @return 最大位数
 */
int MaxBit(int levelQ, int levelP, vector<uint64_t>& qVec, vector<uint64_t>& pVec)
{
    int maxLength = 0;
    for (int i = 0; i < levelQ; i++) {
        uint64_t qi = qVec[i];
        maxLength = max(maxLength, BitsLen64(qi));
    }
    for (int i = 0; i < levelP; i++) {
        uint64_t pi = pVec[i];
        maxLength = max(maxLength, BitsLen64(pi));
    }
    return maxLength;
}

/*
 * @brief ：对64位整数a进行模q减法操作，确保结果在[0, q)范围内
 *
 * @param a：待减法操作的整数
 * @param q：模数
 */
void ComSubQ(uint64_t &a, uint64_t q)
{
    if (a >= q) {
        a -= q;
    }
}

/*
 * @brief ：对64位整数对数组进行模q减法操作，确保结果在[0, q)范围内
 *
 * @param a：待减法操作的数组
 * @param q：模数
 */
void ComSubQNeon(uint64x2_t &a, uint64x2_t q)
{
    uint64x2_t mask = vcgeq_u64(a, q);         // mask = (a >= q) ? 0xFFFFFFFFFFFFFFFF : 0x0
    uint64x2_t sub_value = vandq_u64(mask, q); // mask AND q (preserve q where condition met)
    a = vsubq_u64(a, sub_value);               // a = a - sub_value (only subtract where condition met)
}